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We present a detailed description of the method and results of our calculation of the kaon B parameter using the 
Wilson quark action in quenched QCD at f3 = 5.9 — 6.5. The mixing problem of the As = 2 four-quark operators 
is solved non-perturbatively with full use of chiral Ward identities. We find _B^(NDR, 2GeV) = 0.562(64) in the 
continuum limit, which agrees with the value obtained with the Kogut-Susskind quark action. 



1. Introduction 

The Wilson quark action explicitly breaks chi- 
ral symmetry at finite lattice spacing, which 
causes problems in a number of subjects treated 
by numerical simulations of lattice QCD. For the 
calculation of the kaon B parameter Bk, the 
problem appears as a non-trivial mixing of the 
weak As = 2 four-quark operator of purely left 
handed chirality with those of mixed left-right 
chirality. It has been well known that pertur- 
bation theory does not work effectively for solv- 
ing this mixing problem]!]], and most calcula- 
tions of Bk have tried to resolve the mixing non- 
perturbatively with the aid of chiral perturbation 
theory [0. This method, however, has not been 
successful, since it contains large systematic un- 
certainties from higher order effects of chiral per- 
turbation theory which survive even in the con- 
tinuum limit. 

An essential step for a precise determination 
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of Bk is to control the operator mixing non- 
perturbatively without resort to any effective the- 
ories. The failure of the perturbative approach 
suggests that higher order corrections in terms of 
the coupling constant might be large in the mix- 
ing coefficients. Presence of large 0(a) correc- 
tions in the coefficients is also a possibility. In or- 
der to deal with this problem, the Rome group has 
proposed the method of non-perturbative renor- 
malization(NPR)^], which shows an improve- 
ment of the chiral behavior for the As = 2 op- 
erator [Q. Recently we have proposed an alterna- 
tive non-perturbative method to solve the opera- 
tor mixing problem which is based on the use of 
chiral Ward identities This method fully in- 
corporates the chiral properties of the Wilson ac- 
tion. Our numerical results show that the renor- 
malized As = 2 operator constructed with this 
method has good chiral behavior with the mixing 
coefficients having small momentum scale depen- 
dence. 

The chief findings of our calculation have al- 
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ready been presented in Ref . || . In this report we 
give a detailed description of the implementation 
of our method and the results of our analyses. 

2. Formulation of the method 

Let us consider flavor ST/ (3) chiral variation 
defined by 



8 a ip(x) 
8 a tp(x) 



A a 

tp(x)i— 75, 



(1) 
(2) 



where A a (a = l,---,8) are the flavor matrices 
normalized as Ti(X a X b ) = 2S ab . We consider a set 
of weak operators in the continuum {Oi} which 
closes under flavor chiral rotations 8 a Oi — ic^Oj. 
These operators are given by linear combinations 
of a set of lattice local operators {O a } as Oi = 

^2 a Zi a O a . 

We choose the mixing coefficients Zi a such that 
the Green functions of {O i} with quarks in the 
external states satisfy the chiral Ward identity to 
0(a). This identity can be derived in a standard 
manner || and takes the form given by 



x k 

+ 4(4(0)11^)) 

k 

-i£(6i(0) II i>(Pk)S a i>(Pi)) + 0(a) = 0, 

l kjtl 



(3) 



where pk is the momentum of the external quark, 
Z C J^ and p — (m — 8m) /Z A xt are constants to 
be determined from the Ward identities for the 
axial vector currents (t]], and P a is the pseu- 
doscalar density of flavor a defined by P a (x) = 
^(x)^A°75^(a;). We note that the first term in 
(3) comes from the chiral variation of the Wilson 
quark action and the third represents the chiral 
rotation of the external fields. 

The four-quark operator relevant for E>k 
is given by Ovv+aa = (sj^isj^d) + 
(s7 A i75C?)(s7 / i75(i) where () means color trace. 
To fix the mixing coefficients for the lattice 
four-quark operators, we may choose a par- 
ticular SU(3) flavor chiral rotation to be ap- 
plied for Ovv+aa- Avoiding flavor rotations 



that yield operators which have Penguin contrac- 
tions and hence mix with lower dimension op- 
erators, we employ the A 3 = diag(l,— 1,0) chi- 
ral rotation, under which Oyv+AA and Ova = 
(s7 AI cZ)(s7 At 75<i) form a minimal closed set of the 
operators. 

Since Ovv+aa and Ova are dimension six op- 
erators with As = 2, we can restrict ourselves 
to dimension six operators for the construction of 
the lattice operators corresponding to them. The 
set of lattice bare operators with even parity is 
given by 



VV = (s 7fl d)(s^d) 

AA = (sj f _ l j 5 d)(sj^ 5 d) 

SS = (sd)(sd) 

PP = (s l5 d)(s l5 d) 

TT = (sa^d)(sa^ y d)/2 

and the set with odd parity is 

VA = (s7 M d)(s7 AI 7 5 d) 

SP = (sd)(s^ 5 d) 

TT = (sa f _ t „d)(sa f _ tL ,j 5 d)/2 



(4) 



(5) 



where = [7^,7,^/2. In terms of these opera- 
tors we construct the Fierz eigenbasis, which we 
find convenient when taking fermion contractions 



(G) 



for evaluatin 


g the Green functions in (||), 


On 


= {VV - 


hAA)/2 (+, 


Oi 


= {SS-i 


- TT + PP) 12 (+, 


2 


= {SS- 


- TT/3 + PP) /2 (-, 


0, 


= (VV - 


- AA) /2 + (SS - PP) (-, 


Oi 


= (VV - 


- AA) 12 - (SS - PP) (+, 


5 


= VA 


(+,+) 


O e 


= SP + 


TT/2 (+,-) 


7 


= SP- 


TT/6 (-,-) 



(7) 



Here the first sign after each equation denotes 
the Fierz eigen value and the second the CP^Jll] 
eigen value. We note that the Fierz eigen basis 
we employ is different from the basis chosen by 
the Rome group ||| based on one-loop perturba- 
tion theory. 

The parity odd operators (A3, 7 are CPS odd 
while O5 is CPS even, and hence O5 does not 
mix with Oq 7 under renormalization. Therefore 
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the mixing structure of these operators is given 
by 



OvV+AA 
2Zyv+AA 

Ova 



= O + ziOi 



Z4O4 



'VA 



(8) 



(9) 



where Zyv+AA and Zva are overall renormaliza- 
tion factors. 

Let us consider an external state consisting of 
two s quarks and two d quarks, all having an equal 
momentum p. Under A 3 chiral rotation the Ward 
identity (|§|) for such an external state takes the 
following form: 



Fw+AA = 

-2pZf* (£ x P a (x)\6 vv+AA (0) (ssdd) (p)) 
-(O VA {0)(ssdd)(p)) 

-{±d vv+AA (0)5 3 (ss¥d)(p)) + 0(a) = 0, 

F VA = 

-2pZ%* (E x P a (x)OvA(0) (ssdd) (p) } 

-(^d vv+AA (o)(sfd~d)(p)) 

~(dvA(0)5 3 (S§Td){p)) + 0(a) = 0, 



(10) 



(11) 



where (ssdd)(p) and S 3 (ssdd)(p) repre- 
sent s{p)s(p)d(p)d(p) and s(p)s(p)(d(p)j 5 /2) d(p) 
+s(p)s(p)d(p)(d(p)j 5 /2), respectively. We obtain 
the amputated Green functions for Fyv+AA and 
Fva by truncating the external quark propaga- 
tors according to 



VV+AA,VA 



G- 1 (p)G- 1 ( P )F vv+AAyA G~- 1 (p)G-- 1 (p), 



(12) 



where G" 1 denotes the inverse quark propagator 
with the flavor q. 

Let Pi (i = 0, • • • , 7) be the projection operator 
corresponding to the four-quark operators in the 
Fierz eigenbasis Oi (i = 0, • ■ • , 7). For example, 
we have 

P o^ = £4 hf<" + (7.75) q/3 (7m75) 5A ] ■ (13) 
Since QCD conserves parity one can write 
Tvv+aa/ Zvv+aa — r 5 P 5 , (14) 



v VA /z VA = r p + riPi 



r 4 p 4 . 



(15) 



Expressing Ovv+aa.va in @) in terms of lattice 
operators, we obtain six equations for the five co- 
efficients Zi , • • • , z$ : 

T t = 4+cizH h4z 5 = O(a), i = 0, • • • , 5 (16) 

This gives an overconstrained set of equations, 
and we may choose any five equations to exactly 
vanish to solve for zf. the remaining equation 
should automatically be satisfied to 0(a). We 
choose four equations to be those for i = 1, • • • , 4, 
since 0\, ■ ■ ■ , O4 do not appear in the continuum. 
The choice of the fifth equation, i = or 5, 
is more arbitrary. We have checked that either 
To = or T5 = leads to a consistent result to 
0(a) for zi,---,Zi in the region pa<l. In the 
present analysis we choose r5 = 0. 

The overall factor Zyv+AA is determined by 
the NPR method||] . We calculate the amputated 
Green function. 



G- 1 (p)G 7 1 ( P ) (\6vv+aa(Q) (ssdd) (p)) 
xG g 1 (p)G d 1 ( P ) = f vv+AA (p)P + ■■■ 

and impose the following condition, 

Zvv+aa(p)Z~ 2 (p)T vv+aa (p) = 1, 



(17) 



(18) 



where Z q (p) is the quark wave- function renormal- 
ization factor which is calculated from 



Z q (p) 



^{-iE^mip^G-^p) 
4E„sin 2 (p M ) 



(19) 



We convert the matrix elements on the lattice 
into those of the MS scheme in the continuum 
with naive dimensional regularization (NDR) and 
renormalized at the scale p — 2GeVM: 



B K (NDR,p)=B K (p,l/a 
x "l + ^)f- 4 log^ 

47T \ \p 

where 

(K°\dvv+AA\K°) 



(20) 



14 

y 



;iog2 



B K (p,l/a) 



\(0\A\K°)\< 



(21) 



with p the momentum at which the mixing coeffi- 
cients are evaluated. The axial vector current in 
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Table 1 

Parameters of our simulations. See text for details. 





5.9 


6.1 


6.3 


6.5 


L A x T 


24 a x 64 


32 a x 64 


40 a x 96 


48 a x 96 


#conf. 


300 


100 


50 


24 


skip 


2000 


2000 


5000 


8000 


A 


A 1 COCO 

0.15862 


r\ i c: a no 

U. 15428 


n 1 ci 01 

0.15131 


0.14925 




0.15785 


0.15381 


0.15098 


0.14901 




0.15708 


0.15333 


0.15066 


0.14877 




0.15632 


0.15287 


0.15034 


0.14853 


fitting range (m w ,m p ) 


13-21 


15-25 


18-28 


21 - 31 


fitting range {Bk ) 


19-46 


25-40 


33-64 


36-61 


K c 


0.15986(3) 


0.15502(2) 


0.15182(2) 


0.14946(3) 


a-^GeV) 


1.95(5) 


2.65(11) 


3.41(20) 


4.30(29) 


La(fm) 


2.4 


2.4 


2.3 


2.2 


a Ms( 1 / a ) 


0.1922 


0.1739 


0.1596 


0.1480 


m s a/2 = mdd/2 


0.0294(14) 


0.0198(16) 


0.0144(17) 


0.0107(16) 


S p 2 


1.11 


1.11 


1.15 


1.12 


p* 2 a 2 


0.9595 


0.5012 


0.2988 


0.2056 



the denominator is given by A = Zasj^j^cI with 
Za the renormalization factor determined by the 
NPR method. 

For comparative purpose we also calculate Bk 
with pcrturbative mixing coefficients, for which 
we use the one-loop expression in Ref. || after 
applying a finite correction in conversion to the 
NDR scheme together with the tadpole improve- 
ment with a^jg(l/a). 

Let us remark here that the equations obtained 
in the NPR method corresponds to Ti = 
for i = 1, • • • , 4 in which the contributions of the 
first and the third term in the Ward identity (|3|) 
are dropped. In particular the NPR method ne- 
glects the quark mass contributions coming from 
the first term of (j3|). As may be expected from 
this, the NPR method is equivalent to the Ward 
identities in the limit of large external virtuali- 
ties 1§ . 

3. Parameters of numerical simulation 

Our calculations are made with the Wilson 
quark action and the plaquette action at (3 = 
5.9 - 6.5 in quenched QCD. Table summa- 
rizes our run parameters. Gauge configurations 
are generated with the 5-hit pseudo heat-bath al- 



gorithm. At each value of /3 four values of the 
hopping parameter K are adopted such that the 
physical point for the K meson can be inter- 
polated. The critical hopping parameter K c is 
determined by extrapolating results for ra\ for 
the four hopping parameters linearly in 1/2K 
to m 2 = 0. We take the down and strange 
quarks to be degenerate. The value of half the 
strange quark mass m s a/2 is then estimated from 
m/f/mp = 0.648. 

The physical size of lattice is chosen to be ap- 
proximately constant at La w 2.4fm where the 
lattice spacing is determined from m p = 770MeV. 
To calculate the perturbative renormalization fac- 
tors, we employ the strong coupling constant at 
the scale 1/a in the MS scheme, evaluated by a 
two-loop renormalization group running starting 
from 9^/a) = P av /gf att + 0.0246 with P av the 
averaged value of the plaquette. 

The renormalization factors Zyv+AA, Za and 
the mixing coefficients zi (i = 1, •••,5) are cal- 
culated for a set of external quark momenta 
(i = 1, • • • , ~ 40), which is chosen recursively ac- 
cording to the condition that the (i + l)-th mo- 
mentum p^ l+1 ^a is the minimum number satisfy- 
ing / i+1 ) 2 a 2 > V .pWV with pWa = 2n/T 
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where T denotes the temporal lattice size. We 
employ p* pa 2GeV among pW for the analysis of 
Bk ■ Errors are estimated by the single elimina- 
tion jackknife method for all measured quantities 
throughout this work. 

4. Calculational procedure 

Our calculations are carried out in two steps. 
We first calculate Zyv+AA, Za and z% using the 
quark Green functions having finite space-time 
momenta. For this purpose quark propagators 
are solved in the Landau gauge for the point 
source located at the origin with the periodic 
boundary condition imposed in all four directions. 
For calculating the first term of the Ward iden- 
tity (||), we use the source method JlCfl to insert 
the pseudoscalar density. The point source quark 
propagators are also used for calculating 7r and 
p propagators and to extract their masses from 
them. 

The Bk parameter is extracted from the ratio 
{O K0 {T)6vv + AA(t')O K0 (\)) 



R{t') 



l{O R0 {T)A{t')){A(t')O^ KQ {l)) 



(22) 



— ► B K {p* ,l/a)/L 3 for l«('«r, 

where various operators are defined by Ok ° (t) = 
Ex s(x, t)~f 5 d(x, t), R o (*)=X)a d(x, t)j 5 s(x, t), 
Ovv+AA(t)=J2t=o z w+AAZiO z (x, t) and 
A(t)=J2s Zas(x, t)lH5d(x, t). The contribution 
of the operators Oi (i = 0, • • • , 4) to Bk(p*, 1/a) 
can be obtained similarly from 

{OKo(T)d t (t')Q Ka (l)) 



R\t') 



(23) 



l(0 K o(T)A(f)){A(t')OUl)) 

where Oi(t') = J2s z vv+AAZiOi(x,t'). For cal- 
culation of these ratios we solve quark propaga- 
tors without gauge fixing employing wall sources 
placed at the edges of lattice where the Dirich- 
let boundary condition is imposed in the time 
direction. We obtain Bk at m s /2 by quadrati- 
cally interpolating the data at the four hopping 
parameters. 

5. Results for mixing coefficients 

In Fig. |] we plot a typical result for the mix- 
ing coefficients Z{ (i = 4) as a function 
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Figure 1. Mixing coefficients z\, ■ ■ ■ , za plotted as 
a function of external momentum squared (pa) 2 
for K = 0.15034 at j3 = 6.3. Vertical line corre- 
sponds to p* 2 GeV. 



of the external quark momenta for the case of 
K = 0.15034 at /3 = 6.3. The mixing coeffi- 
cients shows only weak dependence over a wide 
momentum range 0.02<p 2 a 2 <1.0, albeit z\ and 
Z2 have large errors in the small momentum re- 
gion p 2 a 2 <0.1. This enables us to evaluate the 
mixing coefficients with small errors at the scale 
p* pa 2GeV, which always falls within the range 
of a plateau for our runs at (3 = 5.9 — 6.5. 

At the workshop TaleviJllJ] presented a reanal- 
ysis of the results of the Rome group for the mix- 
ing coefficients in the Fierz eigen basis, reporting 
that the momentum dependence in this basis is 
similar to that of our results in Fig. [I]. Their 
simulations are made at /3 = 6.2 with the Clover 
action. For a more detailed comparison, a par- 
allel analysis with the NPR and Ward identity 
methods employing the same quark action on the 
same set of configurations would be desirable. 

In Fig. |^ we compare the mixing coefficients 
Zi (i = 1, • ■ • ,4) evaluated at the scale p* (filled 
symbols) with the perturbative values obtained 
with «j^g(l/a) (open symbols) as a function of 
lattice spacing. We observe that the a depen- 
dence for the mixing coefficients determined from 
the Ward identities is steeper compared to that 
for the perturbative estimates. The magnitude 
of each mixing coefficient for the Ward identity 
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m p a 

Figure 2. Comparison of the mixing coefficients 
Z\ , ■ ■ ■ , Zi evaluated at p* « 2 GeV using the 
Ward identity (WI; solid symbols) and pertur- 
bative (PT; open symbols) methods. The coeffi- 
cients are plotted as a function of m p a. 

method varies nearly in proportion to a, which re- 
duces by 50% between m p a ~ 0.4 and m p a rj 0.2. 

We remark that a large value of Z2 determined 
by the Ward identities sharply contrasts with the 
one-loop perturbative result z 2 = 0. For the other 
coefficients, the perturbative results agree with 
the non-perturbative ones in sign and rough or- 
ders of magnitude. However, they differ in quan- 
titative detail. We find that the magnitudes of 
Z4 are larger than those of z% for all values of j3, 
which is contrary to the perturbative result. 

For our study of the B parameter the mixing 
coefficient Z5 for the parity odd operator Ova is 
not directly relevant. For completeness we plot a 
typical result in Fig. ||. The data shows a scale 
dependence which is stronger than those of Zi (i — 
1, • • • ,4) for parity even operators toward large 
momenta. We do not find this to be particularly 
alarming since Z5 evaluated for a fixed physical 
scale p* approaches unity toward the continuum 
limit as shown in Fig. |[ 

6. Chiral behavior 

Let us examine the chiral property of the op- 
erator Ovv+aa- In Fig. ^ we show the chi- 
ral behavior of the ratio (K °\6 vv+AA \K )/(8/3) 
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Figure 3. Same as Fig. [l] for mixing coefficient 25. 
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Figure 4. Mixing coefficient 2:5 evaluated at p* « 
2 GeV using the Ward identity method as a func- 
tion of m p a. 

/|(0|P|if°)| 2 at (3 = 6.3, where WI stands for 
our method using chiral Ward identities and PT 
for tadpole-improved one-loop perturbation the- 
ory. The solid line represents a quadratic extrap- 
olation of the Ward identity result in the bare 
quark mass m q a = (1/K — l/K c )/2. The extrap- 
olated value at m q = is consistent with zero, 
demonstrating a significant improvement of the 
chiral behavior compared to the perturbative re- 
sult plotted with triangles. 

We plot in Fig. || the values of the ratio extrap- 
olated to m q — as a function of lattice spacing, 
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4- i<oipik°>i 2 




0.00 0.02 
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q c 
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Figure 5. Test of the chiral behavior of 
(K°\d vv+AA \K )/(8/3)/\(0\P\K )\ 2 for the WI 
and PT methods at (3 = 6.3. The solid curve is a 
quadratic extrapolation to the chiral limit. 



where the pseudoscalar density P in the denomi- 
nator is renormalized perturbatively for both WI 
and PT cases (numerical results are given in Ta- 
ble H below). The ratio for the Ward identity 
method becomes consistent with zero at the lat- 
tice spacing m p a<0.3(a<0.08fm). 

In the perturbative approach with the one-loop 
mixing coefficients, chiral breaking effects are ex- 
pected to appear as terms of 0(g 4 ) and 0(a) for 
the Wilson quark action. A roughly linear be- 
havior of our results for the perturbative method 
is consistent with the presence of an O(a) term. 
Making a linear extrapolation to the continuum 
limit a — > 0, we observe that the chiral behavior 
is recovered. This may suggest that the 0(g i ) 
terms in the mixing coefficients left out in the 
one-loop treatment are actually small. 

7. Results for Bk 

We now turn to the calculation of £?^(NDR, 
2GeV). In Fig. [7] we present the ratio R(t') de- 
fined in (^2]) using the mixing coefficients deter- 
mined from the Ward identities for K — 0.15034 
at (3 — 6.3. A good plateau is observed in the 
range 20<i'<70. We make a global fit of the ra- 
tio R(t') to a constant over 33 < t' < 64 for 
this data set. The three horizontal lines denote 
the central value of the fit and a one standard 
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Figure 6. 




'VV+AA 



^°>/(8/3)/|<0|P|iOl 2 



at m q = for the WI and PT methods as a func- 
tion of a. The operators are renormalized at 2 
GeV in the NDR scheme. For both methods we 
use the same P perturbatively corrected with the 
tadpole improvement. The solid line is a linear 
extrapolation to the continuum limit. 

deviation error band. We note that the error of 
the fitted result is roughly equal in magnitude to 
those of the ratio over the fitted range, while we 
would usually expect a smaller error for the fitted 
result. This is because the error of the ratio R(t') 
is governed by those of the mixing coefficients 
(i = l,...,4). 

In Fig. U a representative result for the con- 
tribution of each operator Oi (i = 0,---,4) to 
Bk (p, 1 /a) is plotted as a function of the ex- 
ternal quark momenta, which is obtained by fit- 
ting the ratio R l (t') (i = 0, • • • , 4) of @ with 
a constant over the same fitting range as for 
R(t'). The contributions are nearly independent 
of the external quark momentum in the range 
0.05<p 2 a 2 <1.0 as is expected from the weak scale 
dependence of the mixing coefficients z t shown in 
Fig. A decrease in magnitude of the contribu- 
tions in the small momentum region p 2 a 2 <0.05 
originates from the scale dependence of the over- 
all renormalization factor Zvv+aa/ Z\ estimated 
non-pert urbatively. 

An important observation is that the value of 
Bk (p, f /a) is a result of a large cancellation be- 
tween the amplitudes of ZiOi where each ampli- 
tude of the mixing operator ZiOi (i = 1, • • • , 4) is 
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Figure 7. Ratio R(t') for K = 0.15034 at (3 = 
6.3. Solid lines denote the fitted result and a one 
standard deviation error band. 



comparable to or larger than that of 0$. This is 
the reason why calculations of Bk with the Wil- 
son quark action is difficult. 

We plot the quark mass dependence of 
B K {p*,l/a) at /3 = 6.3 in Fig. |. We observe 
that the results for the PT method seems to di- 
verge toward the chiral limit, while those for the 
WI method stays finite. Interpolating the data 
for the four hopping parameters quadratically to 
m s /2 we obtain the value of Bk(p* , 1/a) at the 
physical point. 

Let us note that the perturbative results have 
quite small errors compared to those of the Ward 
identity method. This is because the mixing co- 
efficients arc definitely given for the perturbative 
method. For a precise determination of Bk with 
the Ward identity method it is of great impor- 
tance to reduce the errors of the mixing coef- 
ficients. For this purpose methods have to be 
devised to effectively compute quark propagators 
for a large set of momenta with good precision. 

Our final results for B K (NDR, 2GeV) obtained 
with (|2(]) are presented in Fig. [l(] as a function of 
lattice spacing. The numerical values are listed 
in Table ||. The method based on the Ward iden- 
tity (WI) gives a value convergent from a lattice 
spacing of m p a ~ 0.3 (Bk ~ 0.6 — 0.8). The large 
error, however, hinders us from making an extrap- 



<K IZ vv+AAZi O,IK> 
-f- l<0IZ A AIK°>| 2 



oO, 

AO, 

oO, 



□□□□□□□□□□□□□□□□□□□□□□ 



0.01 



0.10 



1.00 



10.00 



2 2 

P a 



Figure 8. Contribution of the operators L (i = 
0, • ■ • , 4) to Bk(p, 1/a) with z% determined by the 
WI method for K = 0.15034 at /3 = 6.3. Vertical 
line corresponds to p* w 2 GeV. 

olation to the continuum limit. Since the origin of 
the large error is traced to that of the mixing coef- 
ficients, we develop an alternative method, which 
we refer to as WlySi m which the denominator 
of the ratio for extracting Bk is estimated with 
the vacuum saturation of Ovv+aa constructed 
by the WI method: 



(K°\6 



VV+AA 



^iELo^ l°^ >vs 



(24) 



In this case the fluctuations in the numerator, 
mainly due to those of the mixing coefficients 
Zi, are largely canceled by those in the denom- 
inator, and the resulting error in Bk is sub- 
stantially reduced as apparent in Fig. |l0|. The 
cost is that the denominator of B K S (p* , 1/ a) con- 
tains the contributions of the pseudoscalar den- 
sity (^°|P|0)(0|P|ii: ) besides those of the axial 
vector current (K°\A\0)(0\A\K°), due to which 
the correct chiral behavior of the denominator is 
not respected at a finite lattice spacing. While 
WI and Wlys methods give different results at 
a finite lattice spacing, the discrepancy is ex- 
pected to vanish in the continuum limit. A lin- 
ear extrapolation in a of the Wlys results yields 
Bk(NDR, 2GeV) = 0.562(64), which we take as 
the best value in the present work. This value 
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Table 2 

Results for Bk (NDR, 2GeV) for WI, WI V s and PT methods as a function of f3. Ratio 
(K°\6vv+aa\K )/(8/3)/\(0\P\K q )\ 2 in the chiral limit is also given. 









5.9 


6.1 


6.3 


6.5 


a = 


B K (NBR, 2GeV) 


WI 


+0.38(6) 


+0.68(11) 


+0.69(12) 


+0.72(18) 








Wlys 


+0.168(20) 


+0.288(29) 


+0.342(33) 


+0.360(52) 


+0.562(64) 






PT 


-0.468(14) 


-0.225(22) 


-0.000(21) 


+0.156(40) 


+0.639(76) 


(K"\6vv + AA\K a ) 

!I<o|p|a: >| 2 




WI 

=0 

PT 


-0.0200(39) 
-0.0415(8) 


-0.0068(55) 
-0.0333(10) 


-0.0017(74) 
-0.0240(12) 


+0.006(10) 
-0.0188(17) 


-0.0009(31) 




0.00 



0.02 

m a=(l/K-l/K.)/2 



0.04 



Figure 9. Quark mass dependences of 
B K (p*,l/a) for the WI and PT method at /? = 
6.3. The open symbol is a quadratical interpola- 
tion of the data to m s /2. 

is consistent with a recent JLQCD result with 
the Kogut-Susskind action, (NDR, 2GcV) = 
0.587(7) (17) g§. 

Intriguing in Fig. [lO] is that the perturbative 
calculation (PT), which gives completely "wrong" 
values at a + 1 0, also yields the correct result for 
Bk, when extrapolated to the continuum limit 
a = 0. This is a long extrapolation from negative 
to positive, but the linearly extrapolated value 
Br- (NDR, 2GeV)=0.639(76) is consistent with 
those obtained with the WI or Wlys method. We 
note that a long extrapolation may bring a large 
error in the extrapolated value. 

Finally we mention possible sources of system- 
atic errors in our results from quenching effects 
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Figure 10. Bk (NDR, 2GeV) plotted as a function 
of a for the WI, Wlys and PT methods. The solid 
lines show linear extrapolations to the continuum 
limit. 

and uncertainties for Gribov copies in the Lan- 
dau gauge. With the Kogut-Susskind quark ac- 
tion it has been observed that the error due to 
quenched approximation is small Whether 
this is supported by calculations with Wilson ac- 
tion we must defer to future studies. For the Gri- 
bov problem we only quote an earlier study fllf 
which suggests that ambiguities in the choice of 
the Gribov copies induce only small uncertainties 
comparable to typical statistical errors in current 
numerical simulations. 

8. Conclusions 

Our analysis of Bk demonstrates the effective- 
ness of the method of chiral Ward identities for 
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constructing the As = 2 operator with the correct 
chiral property. We have shown that both Wil- 
son and Kogut-Susskind actions give virtually the 
identical answer for Bk in their continuum limit. 
We may hope that further improvement of our 
simulations, especially the reduction of the errors 
for the mixing coefficients, leads to a precise de- 
termination of Bk with the Wilson quark action. 
The application of this method for calculations of 
Bb is straightforward. 
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